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A general and rigorous methodology to compute the quantum equilibrium isotope effect is de- 
scribed. Unlike standard approaches, ours does not assume separability of rotational and vibrational 
motions and does not make the harmonic approximation for vibrations or rigid rotor approximation 
for the rotations. In particular, zero point energy and anharmonicity effects are described correctly 
quantum mechanically. The approach is based on the thermodynamic integration with respect to 
the mass of isotopes and on the Feynman path integral representation of the partition function. An 
efficient estimator for the derivative of free energy is used whose statistical error is independent of 
the number of imaginary time slices in the path integral, speeding up calculations by a factor of ~ 60 
at 500 K and more at room temperature. We describe the implementation of the methodology in the 
molecular dynamics package Amber 10. The method is tested on three [1,5] sigmatropic hydrogen 
shift reactions. Because of the computational expense, we use ab initio potentials to evaluate the 
equilibrium isotope effects within the harmonic approximation, and then the path integral method 
together with semiempirical potentials to evaluate the anharmonicity corrections. Our calculations 
show that the anharmonicity effects amount up to 30% of the symmetry reduced reaction free energy. 
The numerical results are compared with recent experiments of Doering and coworkers, confirming 
the accuracy of the most recent measurement on 2,4,6,7,9-pentamethyl-5-(5,5- 2 H2)methylene-ll,lla- 
dihydro-12//-naphthacene as well as concerns about compromised accuracy, due to side reactions, 
of another measurement on 2-methyl-10-(10,10- 2 H2)methylenebicyclo[4.4.0]dec-l-ene. 



I. INTRODUCTION 



The equilibrium (thermodynamic) isotope effect (EIE) 
is denned as the effect of isotopic substitution on the equi- 
librium constant. Denoting an isotopolog with a lighter 
(heavier) isotope by a subscript I (h), the EIE is denned 
as the ratio of equilibrium constants 
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where and are molecular partition functions per 
unit volume of reactant and product. We study a specific 
case of EIE - the equilibrium ratio of two isotopomers. 
In this case, the EIE is equal to the equilibrium constant 
of the isotopomerization reaction, 



EIE = K r , 



(1.2) 



where superscripts r and p refer to reactant and product 
isotopomers, respectively. 

Usually, equil ibrium isotope effects are computed only 
approximately.' 1 ' 2 ' 3 ' 4 ' 5 ' '* ' 7 ' 8 ' 9 ' 1 '' ' 11 ' 12 ' " ! In particular, ef- 
fects due to indistinguishability of particles and ro- 
tational and vibrational contributions to the EIE are 
treated separately. Furthermore, the vibrational motion 
is approximated by a simple harmonic oscillator and the 
rotational motion is approximated by a rigid rotor. In 
general, none of the contributions, not even the indistin- 
guishability effects can be separated from the othersJSEH 
However, at room temperature or above the nuclei can 
be accurately treated as distinguishable, and the indis- 
tinguishability effects can be almost exactly described by 



symmetry factors. On the other hand, the effective cou- 
pling between rotations and vibrations, anharmonicity of 
vibrations, and non-rigidity of rotations can in fact be- 
come more important at higher temperatures. For sim- 
plicity, from now on we denote these three effects to- 
gether as "anharmonicity effects" and the approximation 
that neglects them the "harmonic approximation" (HA) . 
In some cases the effects of anharmonicity of the Born- 
Oppenheimer potential surface on the value of EIE can 
be substantial.^ Ishimoto et al. have shown that the 
isotope effect on certain barrier heights^ can even have 
opposite signs, when calculated taking anharmonicity ef- 
fects into account and in the HA.E3 

Our goal is to describe rigorously equilibria at room 
temperature or above. Therefore, two approximations 
that we make are the Born-Oppenheimer approximation 
and the distinguishable particle approximation (we treat 
indistinguishability by appropriate symmetry factors). 
The error due to the Born-Oppenheimer approximation 
was studied for H/D EIE by Bardo and WolfsbergP^ and 
by Kleinman and Wolfsbergf^ and was shown to be of 
the order of 1 % in most studied cases. Since we assume 
that nuclei are point charges, the Born-Oppenheimer ap- 
proximation implies that the potential energy surface is 
the same for the two isotopomers. The differences of 
Born-Oppenheimer surfaces due to differences of nuclear 
volume and qu adrupol e of isotopes can be important for 
heavy element^ 20 ' 21 ' 22 ', but these are not studied in this 
work. 

Symmetry factors themselves result in an EIE equal to 
a rational ratio, which can be computed analytically. In 
order to separate the symmetry contributions from the 
mass contributions to the EIE, it is useful to introduce 
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the reduced reaction free energy, 



II. THE METHODOLOGY 
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where s^ r ' and are the symmetry factors discussed 



in more detail in Sec. Ill To include the effects of quan- 
tization of nuclear degrees of freedom beyond the HA 
rigorously we use the Feynman path integral represen- 
tation (PI) of the partition function. The quantum re- 
duced reaction free energy can then be computed by ther- 
modynamic integration with respect to the mass of the 
isotopes. To compute the derivative of the free energy 
efficiently, we use a generalized virial estimator. The ad- 
vantage of this estimator is that its statistical error does 
not increase with the number of imaginary time slices 
in the discretized path integral. As a consequence, con- 
verged results can be obtained in a significantly shorter 
simulation than with other estimators. 

The ultimate goal would be to combine the path in- 
tegral methodology with ab initio potentials. However, 
since millions of samples are required, the computational 
expense results in the following "compromise:" First, the 
equilibrium isotope effects are computed using ab initio 
potentials, but as usual, within the HA. Then all anhar- 
monicity corrections are computed using the PI method- 
ology, but with semiempirical potentials. In other words, 
we take advantage of the higher accuracy of the ab initio 
potentials to compute the harmonic contribution to the 
EIE and then make an assumption that the anharmonic- 
ity effects are similar for ab initio and semiempirical po- 
tentials. 

After describing theoretical features of the method, we 
apply it to hydrocarbons used in experimental measure- 
ments of isotope effects on [1,5] sigmatropic hydrogen 
transfer reacti ons. Two of them were recently used by 
Doering et aZ p 3 * 24 * who reported equilibrium ratios of 
their isotopomers. This allows us to validate our calcu- 
lations as well as to discuss the apparent discrepancy in 
measurements of Doering et al. from a theoretical point 
of view. 

The outline of the paper is as follows: In Sec. [II] 
we describe a rigorous quantum-mechanical methodol- 
ogy to compute the EIE. Section III presents the [1,5] 



sigmatropic hydrogen shift reactions on which we test 
the methodology, explains how ab initio methods can be 
combined with the PI to compute the EIE, and discusses 
in detail symmetry effects in these reactions. Section |TV| 
explains the implementation of the method in Amber 10 
and describes details of calculations and error analysis of 
our PIMD simulations. Results of calculations are pre- 
sented and compared with experiments in Sec. [V] Section 
VI concludes the paper. 



Thermodynamic integration 

EIE can be calculated by a procedure of thermody- 
namic integration^ with respect to the mass. This 
method takes advantage of the relationship 



EIE = 
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exp 
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where F = — log Q//3 is the (quantum) free energy and 
A is a parameter which provides a smooth transition be- 
tween isotopomers r and p. This can be accomplished, 
e.g., by linear interpolation of masses of all atoms in a 
molecule according to the equation 



(A) = (1 - A) m| r) +Am 
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In contrast to the partition function itself, the integrand 
ofEq. EA\ 
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is a thermodynamic average and therefore can be com- 
puted by either Monte Carlo or molecular dynamics sim- 
ulations. 



Path integral approach 



Classically, the EIE is trivial and Eq. (2.1l can be 



evaluated analytically. When quantum effects are im- 
portant, this simplification is not possible. To describe 
quantum thermodynamic effects rigorously, one can use 
the path integral formulation of quantum mechanics.^ 
In the path integral formalism, thermodynamic prop- 
erties are calculated exploiting the correspondence be- 
tween matrix elements of the Boltzmann o perat or and 
the quantum propagator in imaginary time.LL4I2a j n 
last decades, path integrals proved to be very useful in 
many areas of quantum chemistry, most recently in cal- 
culations of heat capacities,^ rate constants,^ kinetic 
isotope effects,^ or diffusion coefficients.^ 

Let N be the number of atoms, D the number of 
spatial dimensions, and P the number of imaginary time 
slices in the discretized PI (P = 1 gives classical mechan- 
ics, P — > oo gives quantum mechanics). Then the PI 
representation of the partition function Q in the Born- 
Oppenheimer approximation is 
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FIG. 1: Convergence of the generalized virial estimator 
(GVE) and the thermodynamic estimator (TE) as a function 
of the number P of imaginary time slices in the path integral. 
All results obtained by 1 ns long simulations (with the time 
step 0.05 fs) of compound 1-5,5,5-^3 (A = 0) using GAFF 
force field, normal mode PIMD, and Nose-Hoover chains of 
thermostats. 



where r( s ) = nr^r^, . . .,r^J is the set of Carte- 
sian coordinates associated with the s th time slice, and 
$ (jr( s )}) is the effective potential 



FIG. 2: Root mean square errors (RMSEs) of the general- 
ized virial estimator (GVE) and the thermodynamic estima- 
tor (TE) as a function of the number P of imaginary time 
slices in the path integral. All results obtained by 1 ns long 
simulations (with the time step 0.05 fs) of compound 1-5,5,5- 
d 3 (A = 0) using GAFF force field, normal mode PIMD, and 
Nose-Hoover chains of thermostats. Note that the RMSE of 
the GVE is not only non-increasing (as expected from the- 
ory), but in fact decreases slightly with increasing P, which 
is due to the decrease of correlation length. 



Generalized virial estimator 
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The P particles representing each nucleus in P differ- 
ent imaginary time slices are called "beads." Each bead 
interacts via harmonic potential with the two beads rep- 
resenting the same nucleus in adjacent time slices and via 
potential V (r^) attenuated by factor 1/P with beads 
representing other nuclei in the same imaginary time 
slice. By straightforward differentiation of Eq. (1 2. 4| w e 
obtain the so-called thermodynamic estimator (TE), 31 
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A problem with expression (2.7 1 is that its statistical er- 



ror grows with the number of time slices. A similar be- 
havior is a well known property of the thermodynamic 
estimator for energy,^ where the problem is caused by 
the kinetic part of energy. 



The growth of statistical error of the thermodynamic 
estimator for energy is removed by expressing the estima- 
tor only in terms of the potential and its derivatives using 
the virial theorem.^ In our case, a similar improvement 
can be accomplished, if a coordinate transformation, 



XII, (s) 

"V ( r i 



(2.8) 



is introduced into Eq. ( |2.4[ ) prior to performing the 
derivative. Here, the "centroid" coordinate is denned as 



1 p " 1 
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(2.9) 
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In other words, first the "centroid" coordinate r- ; is 
subtracted, and then the coordinates are mass-s caled. 
Resulting generalized virial estimator (GVE) 29 33 takes 
the form 
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Its primary advantage is that the root mean square error 
(RMSE) of the average, er av , is approximately indepen- 
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dent on the number of imaginary time slices P, 

In this equation r s i m denotes the length of the simulation 
and r c is the correlation length. 

The convergence of values and statistical errors of both 
estimators as a function of number P of imaginary time 
slices for systems studied in this paper is discussed in Sec- 
tion [iVj As expected, up to the statistical error they give 
the same values as can be seen in Fig. [T] Nevertheless, 
when quantum effects arc important, and a high value 
of P must be used, the GVE is the preferred estimator 
since it has much smaller statistical error and therefore 
converges much faster than the TE (see Fig. [2]). 

Path integral molecular dynamics 



(PIMC) or path integral molecular dynamics (PIMD). In 
PIMC, gradients of V in Eq. (2.101 result in additional 



calculations since the usual Metropolis Monte Carlo pro- 
cedure for the random walk only requires the values of V. 
This additional cost can be, however, reduced either by 
less frequent sampling, or by using a trick in which the 
total derivative with respect to A (not the gradients!) is 
computed by finite difference! 27 1 29 1 33 1 3 ^ In case of PIMD, 



the presence of gradients of V in Eq. ( 2.10 ) does not slow 



Thermodynamic average in Eq. (2.101 can be eval- 



uated efficiently using the path integral Monte Carlo tious classical momenta p 



down the calculation since forces are already computed 
by a propagation algorithm. Although in principle, a 
PIMC algorithm for a specific problem can always be at 
least as efficient as a PIMD algorithm, in practice it is 
much easier to write a general PIMD algorithm and so 
PIMD is usually the algorithm used in general software 
packages. Since PIMD was implemented in Amber 9,^ 
we have implemented the methodology described above 
for computing EIEs into Amber lO! 3 ^ This implemen- 
tation is what was used in calculations in the following 
sections. In PIMD, equation (2.4 1 is augmented by ficti- 
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The normalization prcfactor C is chosen in such a way 
that the original prefactor C in Eq. (2.4 1 is reproduced 



when the momentum integrals in Eq. (2.121 are evaluated 
analytically. The partition function (2.121 is formally 



equivalent to the partition function of a classical sys- 
tem of cyclic polyatomic molecules with harmonic bonds. 
Each such molecule represents an individual atom in the 
original molecule and interacts with molecules represent- 
ing other atoms via a potential derived from the poten- 
tial of the original molecule.^ This system can be stud- 
ied directly using well developed methods of the classical 
molecular dynamics. 



Low- and high-temperature limits 

It is us eful to see how t he ge neral path integral ex- 
pressions (2.1|, ( 2.4 1 , and (2.101 behave in the low and 



high temperature limits. As temperature decreases, the 
difference of the reduced free energies of isotopomers 
approaches the difference of their zero point energies 
(ZPEs). Therefore, still assuming that indistinguisha- 
bility is correctly described by symmetry factors, the low 
temperature limit of the EIE is equal to 
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where e denotes the ZPE. 

At high temperature, the system approaches its classi- 
cal limit. In this limit, we can set the number of imagi- 
nary time slices P — 1 , and the EIE can be computed an- 
alytically using partition function (2.4), which becomes 
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Again, for isotopomers the mass dependent factors can- 
cel out upon substitution into Eq. (1.2 1, along with 
all other terms except for the integrals. Although the 
Born-Oppcnheimer potential surface is the same for all 
isotopomers, the integrals do not cancel out since the 
integration is not performed over the whole configura- 
tion space, but only over parts attributed to reactants 
or products. Particles are treated as distinguishable in 
the classical limit and the volumes of the configuration 
space corresponding to reactants and products are gener- 
ally different. After the reduction, EIE becomes exactly 
equal to the ratio of the symmetry factors, 



s(p) 



(2.15) 



5 



CD, 



1-5.5.5-rf, 



l-1.5.5-rf 3 



CH^D CHD 



l-l,l-<*, 



1 -J.J-rf, 



1-1.5-,/, 



FIG. 3: Equilibrium of the [1,5] hydrogen shift reaction 
in (3Z)-(5,5,5- 2 H3)penta-l,3-diene (1-5, 5, 5-0(3) and in (3Z)- 
(l,l- 2 H2)penta-l,3-diene (1-1,1-^2). If all contributions ex- 
cept for those due to symmetry factors s a and st were ne- 
glected, one would obtain approximate equilibrium constants 
K\ = 3 and Kt = 2 in both cases. 



III. EQUILIBRIUM ISOTOPE EFFECTS IN 
THREE [1,5] SIGMATROPIC HYDROGEN SHIFT 
REACTIONS 



unwanted side reactions, mainly dimerizations. One 
of our goals is to elucidate this discrepancy from the 
theoretical point of view. 

As can be seen in Figs. [3] and [4] the fi- 
nal equilibrium of the [1,5] hydrogen shift re- 
action of all examined compounds can be de- 
scribed as an outcome of two reactions. The sec- 
ond reaction (leading from deuterio-methyl-dideuterio- 
methylene to dideuterio-methyl-deuterio-methylene in 
tri-deuterated compounds and from dideuterio-methyl- 
mcthylcnc to dcuterio-methyl-deuterio-methylcnc in di- 
deuterated compounds) produces a mixture of species 
that differ by the orientation of the mono-deuterated 
methylene group. Due to a high barrier for rotation of 
this group, the two products cannot be properly sampled 
in a single PIMD simulation. Therefore an additional 
PIMD simulation is necessary, as shown on the example 
of 1-5,5,5- e?3 in Fig. [5] The reduced free energy of the 
second step is then calculated as 
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We examine the EIE for four related compounds. The 
parent compound (3Z)-penta-l,3-diene (compound 1) is 
the simplest molecule to model the [1,5] sigmatropic 
hydrogen shift reaction. Two of its isotopologs, tri- 
deuterated (3Z)-(5,5,5- 2 H 3 )penta-l,3-diene (l-5,5,5-d 3 ) 
and di-deuterated (3Z)-(l,l- 2 H 2 )penta-l,3-diene (1-1,1- 
(I2) (see Fig. [3| were used by Roth and Konig to mea- 
sure an unusually high value of the kinetic isotope effect 
(KIE) on the [1,5] hydrogen shift reaction with respect to 
the substitution of hydrogen by deuterium.^ This result 
pointed to a significant role of tunneling in [1,5] sigma- 
tropic hydrogen transfer reactions. Subsequently, much 
theoretical research was devoted to the study of this re- 
action. Here we calculate final equilibrium ratios of prod- 
ucts of this reaction, which, to our knowledge, were not 
theoretically predicted so far. 

Two other compounds, 2-met hy 1-10- (10,1 0- 
2 H 2 )methylenebicyclo[4.4.0]dec-l-ene (2-l,l-d 2 ) and 
2, 4,6,7, 9-pentamethyl-5-(5,5- 2 H2)methylene-l 1,11a- 
dihydro-12i/-naphthacene (3-1, l-c^) (see Fig. [4]), were 
recently used by Doering et al. to confirm and possibly 
refine the exp erime ntal value of the KIE on the [1,5] 
hydrogen shift J22I23 In contrast to (3Z)-penta-l,3-diene 
(compound 1), where the s-trans conformer incompetent 
of the [1,5] hydrogen shift is the most stable, pentadiene 
moiety in compounds 2 and 3 is locked in the s-cis 
conformation. This not only increases the reaction rate, 
but also rules out the (very small) effect of the EIE on 
the KIE due to the shift in s-cis/s-trans equilibrium. For 
both molecules Doering et al. reported final equilibrium 
ratios of isotopomers. Despite the similarity of com- 
pounds 2 and 3, these ratios are qualitatively different. 
Indeed, one motivation for measuring EIE in compound 
3 was that Doering et al. suspected that in the case 
of 2-1,1- d,2 the equilibrium ratio might be modified by 



where 1/2 is the ratio of symmetry factors and K\ 
and Kf lMB stand for equilibrium constants obtained by 
the second and third PIMD simulations, respectively. To- 
gether they represent the second reaction step. 



Combination of ab initio methods with the PIMD 

Unfortunately, at present the PIMD method cannot 
be used in conjunction with higher level ab initio meth- 
ods, due to a high number of potential energy evaluations 
needed. Semiempirical methods, which can be used in- 
stead, do not achieve comparable accuracy. We therefore 
make the following two assumptions: First, we assume 
that the main contribution to the EIE can be calculated 
in the framework of HA. Second, we assume that selected 
semiempirical methods are accurate enough to reliably 
estimate the anharmonicity correction. The anharmonic- 
ity correction is calculated as 
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With these two assumptions, we can take advantage of 
both PIMD and higher level methods by adding the 
semiempirical anharmonicity correction to the HA result 
calculated by a more accurate method. The HA value 
of AF red is obtained by Boltzmann averaging over all 
possible distinguishable conformations, 



A ^ha = -k B T\n 



s(p) 



where N r is the number of "geometrically different iso- 
mers" of a reactant. By geometrically different iso- 
mers we mean species differing in their geometry, not 
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FIG. 4: Equilibrium of the [1,5] hydrogen shift reaction in 2-methyl-10-(10,10- 2 H 2 )methylenebicyclo[4.4.0]dec-l-ene (2-l,l-d 2 ) 
and in 2,4,6,7,9-pentamethyl-5-(5,5- 2 H2)methylene-ll,lla-dihydro-12i/-naphthacene (3-1, 1-^2) (the locators of positions of 
deuterium atoms in abbreviations are chosen to correspond to locators in (3Z)-penta-l,3-diene). As for compound 1, values of 
equilibrium constants imposed by symmetry are K\ — 3 and K2 = 2. 



species differing only in positions of isotopically substi- 
tuted atoms. Ef is the electronic energy (including nu- 
clear repulsion) of i th isomer, is the symmetry factor 
and <5[; nuc are partition functions of the nuclear motion 



of isotopomers. N p , s^ p \ Qfj nuc denote analogous 
quantities for the product. 



Symmetry factors 

Although symmetry effects can be computed analyti- 
cally, for the reactions studied in this paper they are non- 
trivial and so we discuss them here in more detail. As 
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FIG. 5: Three PIMD simulations used to compute equilib- 
rium ratios of the [1,5] hydrogen shift reaction in (3Z)-(5,5,5- 
2 H3)penta-l,3-diene (1-5,5,5-ds). Half white, half black 
spheres represent deuterium atoms. The methyl group, in 
contrast to methylene, rotates during simulations. 



mentioned above, we are interested in moderate temper- 
atures (above ~ 100 K) where quantum effects might be 
very important but the distinguishable particle approx- 
imation remains valid. In this case, effects of particle 
indistinguishability and of non-distinguishing several, in 
principle distinguishable minima, by an experiment can 
be conveniently unified by the concept of symmetry fac- 
tor. In our setting, we will call "symmetry factor" the 
product 



N 1 



(3.4) 



Here, s cxp refers to the number of distinguishable min- 
ima not distinguished by the experiment and Ui are the 
usual rotational symmetry numbers of symmetric rotors. 
The symmetry numbers are present only if either the 
whole molecule or some of its parts are treated as free 
or hindered classical symmetrical rotors. (In this case, 
the number of minima of the hindered rotor potential is 
not included in s cxp .) The symmetry numbers are not 
present if rotational barriers are so high that the cor- 
responding degrees of freedom should be considered as 
vibrations. 

The concept can be illustrated on an example of the 
mono- and non-deuterated methyl group in a rotational 
potential with three equivalent minima 120° apart. At 
low temperatures, when hindered rotation of the methyl 
group reduces to a vibration, the symmetry factor is de- 
termined only by s cxp . In the case of mono-deuterated 
methyl group there are three in principle distinguishable 
minima corresponding to three rotamers, which are, as 
we suppose, considered to be one species by the observer. 
Therefore s cxp = 3. In the case of non-deuterated methyl 
there is only one distinguishable minimum and s cxp = 1. 
At higher temperatures, when the methyl group can be 
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treated as a hindered rotor, its contribution to the sym- 
metry factor is determined only by rotational symmetry 
number cr, where a = 1 for a mono-deuterated methyl 
group and a = 3 for a non-deuterated methyl group. 
From definition (3.4 1, it is clear that the high and low 



temperature pictures are consistent and give the same 
ratios of symmetry factors. 

Partition functions Q^ r > and needed in calculation 
of the reduced free energy in Eq. ( |1.3| are computed as 
sums of partition functions of sW' ox p and s(p)> cx p iso- 
topomers. 



IV. COMPUTATIONAL DETAILS 



Ab initio, density functional, and semiempirical 
methods 



In this subsection, we will discuss the accuracy of four 
electronic structure methods used in our calculations . Ab 
initio MP 2 and the B98 density functional methodJ^ESl 
both in combination with the 6-311+(2df,p) basis set, 
were used for calculations within the HA. Semiempirical 
AMIES and SCC-DFTB^ methods were used in both HA 
and PIMD calculations. This allowed us to compute the 
error introduced by the HA. 

Aside from symmetry factors, the EIE is dominated 
by vibrational contributions. Therefore we concentrate 
mainly on the accuracy of HA vibrational frequencies. 
According to Merrick et aZ.pU who tested the perfor- 
mance of MP2 and B98 methods by means of compar- 
ison with experimental data for a set of 39 molecules, 
RMSE of ZPEs is 0.46 kJ • mol" 1 at MP2/6-311+(2df,p) 
and 0.31 kJ -mor 1 at B98/6-311+(2df,p) level of the- 
ory. Appropriate ZPE scaling factors are equal to 0.9777 
and 0.9886 respectively. Corresponding RMSEs of fre- 
quencies are 40 cm -1 and 31 cm -1 . Therefore, a slightly 
higher accuracy can be expected from the B98 functional. 
The accuracy of both semiempirical methods is signif- 
icantly worse than the accuracy of higher level ab ini- 
tio methods. According to Witek and Morokuma, the 
RMSE of AMI frequencies in comparison to experimental 
values for 66 molecules is 95 cm' 1 with frequency scaling 
factor equal to 0.9566M^The error of vibrational frequen- 
cies obtained using SCC-DFTB depends on parametriza- 
tion. We tested two parameter sets: the original SCC- 
DFTB parametrizatiorP^ and the parameter set opti- 
mized with respect to frequencies by Malolepsza et aZp^ 
(further designated SCC-DFTB-MWM). The error of vi- 
brational frequencies calculated with the original SCC- 
DFTB parameters was studied by Kriiger et aZP^ The 
mean absolute deviation from the reference values calcu- 
lated at BLYP/cc-pVTZ for a set of 22 molecules was 
75 cm -1 . The error of the reference method itself, as 
compared to an experiment with a slightly smaller set of 
molecules was 31 cm . In the above mentioned study 
performed by Witek and Morokuma, RMSE of 82 cm -1 
with scaling factor 0.9933 was obtainedPHl The SCC- 



DFTB-MWM mean absolute deviation of experimental 
and calculated frequencies for a set of 14 hydrocarbons 
is indeed better and is equal to 33 cm -1 instead of 59 
cm -1 for the original parametrizationP^l 

The suitability of AMI, SCC-DFTB, SCC-DFTB- 
MWM, and several other semiempirical methods for our 
systems was tested by comparison of EIE values in the 
HA and of potential energy scans with the correspond- 
ing quantities computed with MP2 and B98. Results 
of this comparison are presented in the appendix. The 
AMI method is shown to reproduce ab initio EIEs in the 
HA very well, but fails to reproduce scans of potential 
energies of methyl and vinyl group rotations. There- 
fore, in addition to AMI method we used the SCC- 
DFTB method, which is, among the semiempirical meth- 
ods tested by us, the best in reproducing ab initio poten- 
tial energy scans. On the other hand, compared to AMI, 
SCC-DFTB gives a worse EIE in the HA. 



Statistical errors, convergence, and parameters of 
PIMD simulations 



Statistical RMSEs of averages of PIMD simulations 
were calculated from the equation 



(T. 



(4.1) 



where a is the RMSE of one sample. Correlation lengths 
t c were estimated by the method of block averages.^ At 
constant temperature, the correlation length decreases 
with increasing number of imaginary time slices. Also, 
for constant number of imaginary time slices, the cor- 
relation length decreases with increasing temperature. 
Finally, the correlation length stays approximately con- 
stant at different temperatures, if the number of imag- 
inary time slices is chosen so that AF is converged to 
approximately the same precision. In our systems, cor- 
relation length is close to 3.5 ps. 

The EIE was studied at four different temperatures: 
200, 441.05, 478.45, and 1000 K using the normal mode 
version of PIMD.HSl To control temperature, the Nose- 
Hoover chains with four thermostats coupled to each PI 
degree of freedom were used.^1. 

Different numbers P of imaginary time slices had to 
be used at different temperatures, since at lower tem- 
peratures quantum effects become more important, and 
the number of imaginary time slices necessary to main- 
tain the desired accuracy increases. To examine the re- 
quired value of P as a function of temperature we used 
the GAFF force fieldPS Whereas the accuracy of vibra- 
tional frequencies calculated using the GAFF force field 
is relatively low [RMS difference of GAFF and B98/6- 
311+(2df,p) frequencies of compound 1 is equal to 125 
cm -1 ], the potential should be realistic enough for the 
assessment of the convergence with respect to P. For ex- 
ample, the difference of potential energies AE between 
s-cis and s-trans conformations of compound 1 is 4.3 
kcal • mol -1 as compared to 2.7 kcal • mol -1 obtained 
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by MP2/6-311+(2df,p) or 3.5 kcal ■ mor 1 obtained by 
B98/6-311+(2df,p). 

To check the convergence at 478.45 K, we calculated 
values of the integral in Eq. (2.1l for the deuterium trans- 



The integral in Eq. (2.1l was calculated using the 



fer reaction in l-5,5,5-d 3 with 40 and 48 imaginary time 
slices (using Simpson's rule with 5 points). Their differ- 
ence is equal to 0.00005 ± 0.00040 kcal ■ mor 1 . This is 
less than the statistical error of the calculation on the 
model system, which itself is smaller than the error of 
production calculations, since the model calculation was 
10 times longer than the longest production calculation. 
Therefore, taking into account the accuracy of produc- 
tion calculations, P = 40 can be considered the con- 
verged number of imaginary time slices. This is further 
supported by the observation of the convergence of the 
single value of the GVE (pUOj at A = 0. The relative 
difference of GVE values obtained with 40 and 48 imag- 
inary time slices is equal to 0.22 ± 0.02 %, whereas the 
difference between 40 and 72 imaginary time slices equals 
to 0.39 ± 0.04 %. The discretization error is asymptot- 
ically proportional to P~ 2 . By fitting this dependence 
to the calculated values, we estimated the difference be- 
tween the values for P = 40 and for the limit P — > oo to 
be less than 0.6 %. This is only three times more than 
the difference between P = 40 and P = 48. Therefore, if 
the integral converges similarly as the GVE, we can use 
the aforementioned difference between 40 and 48 imag- 
inary time slices as the criterion of convergence. The 
convergence of the GVE with the number of imaginary 
time slices is displayed in Fig. [I] Since 441.05 K is close 
enough to 478.45 K we used the same value of P at this 
temperature. To check the convergence at 200 and 1000 
K, we observed only the convergence of the single value 
of the GVE. At 200 K the relative difference of GVE val- 
ues at A = obtained with P = 72 and P = 80 is 0.03 
± 0.07 %. Based on the comparison with the previous 
result, P = 72 is considered sufficient. At 1000 K, the 
relative difference of GVE values at A = for P = 24 
and P = 32 equals to 0.1 ± 0.02 %, so that 24 imaginary 
time slices are used further. 

The time step at 441.05 and 478.45 K was 0.05 fs to 
satisfy the requirement of energy conservation. At 200 
and 1000 K, a shorter step of 0.025 fs was used due to 
the increased stiffness of the harmonic bonds between 
beads at 200 K and due to the increased average kinetic 
energy at 1000 K. The simulation lengths differed for 
different molecules. For both isotopologs of compound 1 
simulation length of 1 ns ensured that the system prop- 
erly explored both the s-trans and the s-cis conforma- 
tions. Convergence was checked by monitoring running 
averages and by comparing the ratio of the s-trans and s- 
cis conformers with the ratio calculated in the HA. The 
length of converged PIMD simulations of compound 2 
was 500 ps. Convergence was checked again using run- 
ning averages and by visual analysis of trajectories to 
ensure that the system properly explored all local min- 
ima. For compound 3, the simulation length was 400 
ps. 



Simpson's rule. Using the AMI potential, the GVE 
was evaluated for five values of A, namely for A = 
0.0, 0.25, 0.5, 0.75, 1.0. Convergence was checked by com- 
parison with values obtained using the trapezoidal rule. 
Since the dependence of the estimator on the parameter 
A is almost linear, the difference between the two results 
remained well under the statistical error. Using the SCC- 
DFTB potential, the dependence on A was less smooth. 
As a result, nine equidistant values of A were needed to 
achieve similar convergence. In one case, as many as 17 
values of A had to be used. 



Amber 10 implementation 



The PIMD calculations were performed using Amber 
lOP^ The part of Amber 10 code, which computes the 
derivative dF (A) /d\ with respect to the mass was im- 
plemented by one of us and can be invoked by setting 
ITIMASS variable in the input file. Several posible ways 
to compute the derivative are obtained by combining 
one of two implementations of PIMD in sander (either 
the multisander implementation or the LES implementa- 
tion) with either the Nose-Hoover chains of thermostats 
or the Langevin thermostat, with the normal mode or 
"primitive" PIMD, and with the TE or GVE. Calculat- 
ing the value of dF (A) /d\ for compound 1-5, 5, 5-^3 us- 
ing GAFF force field, at A = 0, T = 478.45 K, and for 
several values of P, we confirmed that all twelve possible 
combinations give the same result. For example, Fig. [I] 
shows the agreement of the GVE and TE. Nevertheless, 
the twelve methods differ by RMSEs of dF (A) /dX (due 
to different statistical errors of estimators and correla- 
tion lengths) and by computational costs. As expected, 
the most important at higher values of P is the differ- 
ence of statistical errors of GVE and TE. Figure [2] com- 
pares the dependence of the RMSE of the GVE and TE 
on the number P of imaginary time slices. For P = 64, 
the converged simulation with GVE is approximately 100 
times faster than with TE. Less significant differences of 
RMSEs are due to differences in correlation lengths. As 
expected, for primitive PIMD with Langevin thermostat 
the correlation length depends strongly on collision fre- 
quency 7 of the thermostat. The correlation length is 
approximately 450-500 fs for 7 = 0.3 ps" 1 , falling down 
quickly to 120-150 fs for 7 = 3 ps" 1 and to 5 fs for 
7 = 300 ps" 1 , then rising slowly again. This can be com- 
pared to correlation length of 10-20 fs of the primitive 
PIMD thermostated by Nose-Hoover chains of 4 ther- 
mostats per degree of freedom. A smaller difference can 
be found between correlation lengths of the normal mode 
(3-8 fs) and primitive PIMD (10-20 fs). 
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1 st step AF rcd AAF anharm 2 nd step AF rcd AAF anharm 
(3Z)-(5,5,5- 2 H 3 )penta-l,3-diene (l-5,5,5-d 3 ) 

AMI (PIMD) 0.0395 -0.0041±0.0009 -0.0154 0.0022±0.0007 

SCC-DFTB (PIMD) 0.1245 -0.0039±0.0007 -0.0616 0.0026±0.0005 

B98 (HA) + AAF£Mi rm 0.0587 -0.0248 

MP2 (HA) + AAgggy™ 0.0770 -0.0338 

(3Z)-(l,l- 2 H 2 )penta-l,3-diene (l-l,l-d 2 ) 

AMI (PIMD) -0.0283 0.0063±0.0009 0.0191 -0.0023±0.0007 
SCC-DFTB (PIMD) -0.1142 0.0049±0.0006 0.0610 -0.0017±0.0005 
B98 (HA) + AAFXSii rm -0.0466 0.0282 
MP2 (HA) + AAFamT™ -0.0645 0.0372 



TABLE I: Reduced free energies AF rcd [kcal • moP 1 ] and anharmonicity corrections AAF anharm [kcal ■ moP 1 ] of [1,5] hydrogen 
shift reactions in 1-5,5, 5-d 3 and in l-l,l-d 2 at 478.45 K. For AMI and SCC-DFTB methods, values calculated by PIMD 
are listed followed by the anharmonicity correction obtained as the difference of the PIMD and HA value. For B98 and MP2 
methods, HA values corrected by the AMI anharmonicity correction are listed. The first reaction step in the case of tri- 
deuterated compound leads from 1-5,5,5-^3 to 1-1,1,5-rfs, which is also the reactant of the second reaction step leading to 
1-1, 5, 5-^3. In case of di-deuterated compound the sequence is: l-l,l-d 2 , l-5,5-d 2 and l-l,5-d 2 . 



Used software 

All PIMD calculations were performed in Amber 10P^ 
All MP2 and B98 calculations as well as AMI and PM3 
scmicmpirical calculations in the HA were done in Gaus- 
sian 03 revision EOlP DFTB and SCC-DFTB calcula- 
tions in the HA used the DFTB+ code, version 1.0. l.^ 
SCC-DFTB harmonic frequencies were computed numer- 
ically using analytical gradients provided by the DFTB+ 
code. The step size for numerical differentiation was set 
equal to 0.01 A. This value was also used by Kriiger et 
aPsl in their study validating the SCC-DFTB method 
and their frequencies differed from purely analytical fre- 
quencies of Witek et alW by at most 10 cm . To di- 
agonalize the resulting numerical Hessian, we used the 
f ormchk utility included in the Gaussian program pack- 
age. 



V. RESULTS 

A. (3Z)-penta-l,3-diene 

(3Z)-pcnta-l,3-diene (compound 1) is the simplest of 
examined molecules. Its non-deuterated isotopolog has 
three distinguishable minima: the s-trans conformer, 
which is the global minimum and two s-cis conformers re- 
lated by mirror symmetry. Strictly speaking, in pentadi- 
ene the s-cis species have actually gauche conformations 
due to sterical constraints. In their original experiments, 
Roth and Konig studied two isotopologs, 1-5,5,5- e?3 and 
l-l,l-d 2 - The EIEs of both isotopologs were computed 
using the PIMD methodology of Sec. [rfjat 478.45 K. Re- 
sulting reduced free energies are listed in Table [T] Note 
that the anharmonicity correction is very similar for AMI 
and SCC-DFTB methods, even though the main part of 
AF rcd obtained in the HA (AFf^) substantially differs 



for the two methods. This indicates that the corrections 
are fairly reliable in this case and can be used to cor- 
rect results of higher level methods obtained in the HA. 
The anharmonicity correction is as large as 20 % of the 
final value of the reduced free energy. Unfortunately, 
the difference between MP2 and B98 in the HA is still 
approximately four times larger than the anharmonicity 
correction. 



1-5,5, 5-d 3 l-l,l,5-d 3 l-l,5,5-d 3 



AMI (PIMD) 


0.103 


0.296 


0.601 


SCC-DFTB (PIMD) 


0.108 


0.285 


0606 


B98 (HA) + AAFlST rm 


0.104 


0.293 


0.602 


MP 2 (HA) + AAFJ'Sf™ 


0.105 


0.291 


0.604 



l-l,l-d 2 l-5,5-d 2 l-l,5-d 2 



AMI (PIMD) 


0.099 


0.305 


0.597 


SCC-DFTB (PIMD) 


0.093 


0.315 


0.591 


B98 (HA) + AAFlST rm 


0.097 


0.307 


0.596 


MP 2 (HA) + AAFXmi™ 


0.096 


0.309 


0.595 



TABLE II: Equilibrium ratios of [1,5] hydrogen shift reactions 
of l-5,5,5-d 3 and 1-1,1-da at 478.45 K. 



Reduced free energies suggest a general preference 
valid for both tri- and di-deuterio species: Namely, deu- 
terium, compared to hydrogen, prefers sp 3 carbon of the 
methyl group to sp 2 carbon o f the vin yl group. This was 
also observed experimentally. ^ 54 * 55 ! As can be seen in 
the table, the preference is present already in the HA. 
At first sight this preference can be counter-intuitive: If 
we confine ourselves only to the most energetical C-H 
(or C-D) bond stretching modes and suppose that force 
constants do not change significantly upon substitution 
with deuterium, the deuterium should prefer the stiffcst 
bonds. This is because more energy can be gained by 
substituting the stiffer sp 2 C-H bonds, assuming approx- 
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1 st step AF rcd AAF anlmrm 2 nd step AF rcd AAF anharm 

AMI (PIMD) -0.0337 0.0102±0.0013 0.0230 -0.0036±0.0010 

B98 (HA) + AAFamT"" -0.0496 0.0302 

MP2 (HA) + AAJami™ -0.0674 0.0392 

TABLE III: Reduced free energies AF Icd [kcal ■ moR 1 ] and anharmonicity corrections AA.F anharm [kcal ■ moR 1 ] of the [1,5] hy- 
drogen shift reaction in 2-methyl-10-(10,10- 2 H2)methylenebicyclo[4.4.0]dec-l-ene (compound 2) at 441.05 K. For AMI method, 
values calculated by PIMD are listed followed by the anharmonicity correction obtained as the difference of PIMD and HA 
values. For B98 and MP2 methods, only HA values corrected by the AMI anharmonicity correction are listed. The first 
reaction step leads from 2-l,l-d 2 to 2-5,5-d 2 , which is also the reactant of the second reaction step leading to 2-l,5-d 2 . 



imately the same change in reduced mass after substi- 
tution. (Recall that the energy of a vibrational mode is 
proportional to y/kjji where k is the force constant and 
\i the reduced mass.) Considering the stretching modes 
only, this would be the case for isotopologs of compound 
1. Nevertheless, taking into account also the bending and 
torsional vibrational modes, gains on the sp 3 C-H bond 
side will dominate. (This "counter-intuitive" preference 
of heavier isotope in "softer" bonds is quite common. For 
examples see the inverse H/D equilibrium isotope effect 
in oxidati ve addition reactions of H 2 to transition metal 
complexes) 56 1 57 1 58 1 5 ^ or the inverse 16 0/ ls O isotope ef- 
fect in metal mediated oxygen activation reactionP21) As 
already stressed, final equilibrium concentrations are de- 
termined mainly by the symmetry factors and the afore- 
mentioned deuterium sp 3 to sp 2 preference manifests it- 
self only in a small modification of the symmetry deter- 
mined rational ratios, as seen in Table [TTJ. 
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FIG. 6: Reduced reaction free energies of the first step of [1,5] 
hydrogen shift reaction in (3Z)-(5,5,5- 2 H3)penta-l,3-diene (1- 
5,5,5-d 3 ) calculated in the HA as the Boltzmann average of 
all s-trans and s-cis isomers. 



Temperature dependence of the reduced free energy 

Temperature dependence of the reduced free energy 
for the first reaction step of hydrogen shift in 1-5, 5, 5-^3 
is depicted in Fig. [6] Analogous temperature depen- 
dence for all other studied reactions is very similar to 
the dependence in this figure. At very low temperature, 
the reduced reaction free energy approaches the differ- 
ence of ZPEs in accordance with Eq. (2.13). Absolute 



value of the reduced free energy is maximal at temper- 
atures around 200 K. At temperatures around 400-500 
K, where most measurements took place, the value of 
AF red is (by chance) close to the difference of ZPEs. 
At high temperatures, AF red goes to zero in accordance 
with the high temperature limit (2.15) discussed above, 
which is valid also in the HA as can be shown using 
the Teller- Redlich theorem ! 60 * 61 ! The temperature depen- 
dence of the anharmonicity correction was examined us- 
ing the AMI semiempirical method at temperatures 200 
K, 478.45 K and 1000 K. The correction for the first 
reaction step of 1-5,5,5- d 3 is equal to -0.0042 ± 0.0009, 
-0.0041 ± 0.0009 and -0.0035 ± 0.0008 kcal-moR 1 , 
respectively. Therefore, taking into account statistical 
errors, the correction stays approximately constant over 
the wide temperature range. Since the value of AF rcd is 
decreasing in the region from 200 to 1000 K, the relative 



importance of the correction is increasing. 



B. 2-methyl-10-methylenebicyclo[4.4.0]dec-l-ene 



2-l,l-d 2 2-5,5-d 2 2-l,5-d 2 



AMI (PIMD) 


0.098 


0.306 


0.595 


B98 (HA) + AAi^jjJf™ 


0.097 


0.308 


0.595 


MP2 (HA) + AAFl^ h r m 


0.096 


0.310 


0.594 


experimental (series 1) 


0.108 


0.328 


0.564 


experimental (series 2) 


0.114 


0.314 


0.572 



TABLE IV: Equilibrium ratio of the [1,5] hydrogen shift reac- 
tion in dideuterated compound 2 at 441.05 K. Experimental 
series 1 and 2 were obtained by two different methods of anal- 
ysis of the NMR spectrumP3 



Compound 2-l,l-d 2 was used relatively recently by 
Doering and ZhafP^ to refine the original results of Roth 
and Konig. Doering and Zhao reported the equilibrium 
concentrations at three temperatures, from which we 
have chosen the lowest, T = 441.05 K. Since AMI and 
SCC-DFTB methods gave similar anharmonicity correc- 
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1 st step AF rcd AAF anharm 2 nd step AF Icd AAF anharm 

AMI (PIMD) -0.0439 0.0160±0.0015 0.0265 -0.0080±0.0011 
B98 (HA) + AAFjSi h r m -0.0673 0.0376 

TABLE V: Reduced free energies AF rcd [kcal ■ mo-P 1 ] and anharmonicity corrections AAF anharm [kcal ■ moP 1 ] of the [1,5] hy- 
drogen shift reaction in 2,4,6,7,9-pentamethyl-5-(5,5- 2 H2)methylene-ll,lla-dihydro-12//-naphthacene (compound 3) at 441.05 
K. For AMI method, values calculated by PIMD are listed followed by the anharmonicity correction obtained as the difference 
of PIMD and HA values. For B98, HA values corrected by the AMI anharmonicity corrections are listed. The first reaction 
step leads from 2-1, 1-^2 to 2-5,5-^2, which is also the reactant of the second reaction step leading to 2-1,5-^2. 



tions for compound 1, we used only the AMI method 
in this case. Four minima were found. Energy differ- 
ences between the global minimum and local minima at 
the B98/6-311+(2df,p) level of theory are 4.8, 7.5 and 
9.5 kcal-mol -1 . Resulting AF rcd calculated according 
to Eq. (3.3 1 are listed in Table III As expected, they 
are similar to AF rcd of 1-1,1-^2- Anharmonicity correc- 
tions changed more and they are about 60 % higher in 
the absolute value. 

Calculated equilibrium ratios are listed in Table |IV] to- 
gether with experimental ratios reported by Doering and 
Zhao. Theoretical and experimental ratios differ sub- 
stantially, which suggests that side reactions suspected 
by Doering and Zhao had indeed occurred and influenced 
the accuracy of results of their study. 



C. 2, 4,6,7, 9-pentamethyl-5-methylene- 11,1 la- 
dihydro-12if-naphthacene 

In order to suppress unwanted side reactions suspected 
for compound 2, 2 ? Doering and Keliher further developed 
the model compound into 3-1,1- d 2 by adding methyl 
substituted aromatic rings on both sides of the cyclic 
partJ^S With this compound they obtained the same equi- 
librium ratios for all temperatures they had examined.^ 
Because of this and because the temperature 441.05 K we 
have chosen for our analysis of compound 2 differs from 
one of temperatures used in Ref. HUby less than 2 K, we 
have decided to use this temperature also for compound 
3. In the HA, only the B98 density functional method 
was used, due to a considerable size of the molecule. Be- 
cause of the increased rigidity imposed by aromatic rings 
on the sides of the original bi-cyclic compound, only three 
distinct minima were found. (We neglected several pos- 
sible orientations of two methyl groups distant from the 
reaction site, which hardly affect the final result.) At the 
B98/6-311+G(2df,p) level, the local minima have ener- 
gies 1.3 and 5.6 kcal • mol -1 above the global minimum. 
The AMI method gives the opposite order of the first and 
second lowest minima. Reduced free energies and anhar- 
monicity corrections obtained using the AMI method are 
listed in Table ED 

Values of both AF red and the anharmonicity correc- 
tions arc again qualitatively similar (but higher in ab- 
solute values) to those of the smaller and less strained 
compound 2. Here, values of anharmonicity corrections 





3-l,l-d 2 


3-5,5-d2 


3-l,5-d 2 


AMI (PIMD) 


0.098 


0.308 


0.595 


B98 (HA) 


0.095 


0.312 


0.593 


B98 (HA) + AAFXmi™ 


0.096 


0.310 


0.594 


experimental" 


0.098 


0.308 


0.594 



"Ref. 1231 

TABLE VI: Equilibrium ratios of the [1,5] hydrogen shift re- 
action in compound 3 at 441.05 K. 



reached approximately 30 % of the values of reduced free 
energies of both reaction steps. Resulting equilibrium ra- 
tios together with their experimental values can be seen 
in Table |VI] Agreement of the theoretical prediction with 
the experimental result is very good. An uncorrected B98 

to demonstrate 



HA value is also included in Table VI 
how the anharmonicity correction modifies the HA equi- 
librium ratio. Surprisingly, the direct AMI PIMD calcu- 
lation is closest to the experimental value, but this ought 
to be ascribed to a fortunate coincidence, considering the 
aforementioned accuracy of electronic structure methods 
used in our study. 



VI. DISCUSSION AND CONCLUSIONS 

To conclude, the combination of higher level methods 
in the HA with PIMD using semiempirical methods for 
the rigorous treatment of effects beyond the HA proved 
to be a viable method for accurate calculations of EIEs. 
Using the generalized virial estimator for the derivative of 
the free energy with respect to the mass we were able to 
obtain accurate results at lower temperatures in reason- 
able time (~ 60 times faster than with thermodynamic 
estimator), since the statistical error is independent on 
the number of imaginary time slices. Two semiempirical 
methods, AMI and SCC-DFTB, were used for calcula- 
tion of the anharmonicity correction, both giving very 
similar results. Calculations showed, that the anhar- 
monicity effects account up to 30 % of the final value 
of the reduced free energy of considered reactions. The 
anharmonicity correction always decreases the absolute 
value of the reduced reaction free energy. This is consis- 
tent with the qualitative picture in which higher vibra- 
tional frequencies of hydrogens are lowered by the an- 
harmonicity of a potential more than lower frequencies 
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of deuteriums. This in turn is due to a higher amplitude 
of vibrations of lighter hydrogen atoms. The lower dif- 
ference between frequencies of unsubstituted and deuter- 
ated species results in the lower absolute value of the 
reduced reaction free energy. Unfortunately, the inaccu- 
racy of the ab initio electronic structure methods used 
in our study is still of at least the same order as the 
anharmonicity corrections. 

For isotopologs of compound 1, we predicted equilib- 
rium ratios and free energies of the [1,5] sigmatropic hy- 
drogen shift reaction. A comparison with experimental 
results was not possible due to a low precision of the orig- 
inal measurement. For compound 2, the disagreement 
between theoretical and experimental data supports the 
suspicion by authors of the measurement that the ac- 
curacy of their results was compromised by dimerization 
side reactions. On the other hand, the agreement of theo- 
retically calculated ratios with experimental observations 
in the case of compound 3 suggests that the isolation 
of the [1,5] hydrogen shift reaction from disturbing in- 
fluences was successfully achieved and the observed EIE 
and KIE can be considered reliable. 
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APPENDIX A: EXAMINATION OF 
SEMIEMPIRICAL METHODS USED IN PIMD 
SIMULATIONS 

To determine the suitability of semicmpirical methods 
used in our PIMD calculations, we first compared val- 
ues of the EIE in the HA. The B98 and MP 2 methods 
served as a reference. As can be seen from Fig. 6 in 
the paper, which shows the temperature dependence of 
AF led of the first step of deuterium transfer reaction 
in l-5,5,5-ei3, the difference between B98 and MP2 at 
temperatures below 500 K is close to 0.02 kcal/mol. So 
is the difference between AMI and B98 methods. On 
the other hand, the SCC-DFTB method clearly overes- 
timates the extent of EIE compared to both higher level 
methods. A very similar trend was observed in all exam- 
ined reactions. Other semiempirical methods tested were 
PM3^1 and SCC-DFTB-MWM, which are not included 
in Fig. 6 in the paper for clarity. The PM3 method over- 
estimates the EIE similarly to SCC-DFTB, whereas the 
SCC-DFTB-MWM curve is somewhat closer to ab initio 
curves than the SCC-DFTB one. To conclude, from this 
point of view AMI is the preferred semiempirical method. 

During simulations, the pentadiene molecule often 
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FIG. 7: Relaxed potential energy scan of the methyl rota- 
tion in s-trans (3Z)-penta-l,3-diene (compound 1). For MP2, 
B98, DFTB, SCC-DFTB and SCC-DFTB-MWM only posi- 
tions and potential energies of minima and maxima are indi- 
cated. 



passes two potential energy barriers. These are the bar- 
rier for the hindered rotation of the methyl group and 
the barrier for the rotation of the vinyl group, which con- 
nects s-trans and s-cis conformations. Relaxed potential 
energy scans of these two motions were employed as the 
second criterion to assess the relevancy of semiempirical 
methods. Methods tested were MP2, B98^ AMI, SCC- 
DFTB, SCC-DFTB-MWM, PM3, RMlp^J PM3CARB- 
lilPDDG/PM3,IS]and PM6P Potential surface scans 
with PM3CARB1, RM1, PM3/PDDG methods were 
calculated using the public domain code MOPAC 6. 
PMG potential surface scans were performed in MOPAC 
2007p3 Results for the methyl group rotation are shown 
in Fig. [7] The height of the AMI barrier is only 0.005 
kcal • mol . Moreover, positions of minima do not agree 
with B98 and MP2. On the other hand, the SCC-DFTB 
method matches higher level methods closely. From other 
semiempirical methods PM3 performs best in this aspect. 
The height of the barrier is relatively well reproduced also 
by PDDG/PM3, RM1 and SCC-DFTB-MWM methods, 
but positions of extrema of the potential energy surface 
are incorrect. Figure [8] shows potential energy scans of 
the s-trans /s-cis rotation of the vinyl group. Again, 
the SCC-DFTB method matches higher level methods 
closely. All other methods (with the exception of DFTB) 
give too low barrier heights as well as too low energy dif- 
ferences between s-trans and s-cis conformations. Also 
note that the potential energy surfaces of PM3 and re- 
lated methods (PDDG/PM3 and PM3CARB-1) are not 
smooth in the gauche region. This peculiarity of the PM3 
potential surface can be seen also in the potential surface 
scan performed by Liu et alW^ Based on these results we 
concluded that none of the semiempirical methods except 
for SCC-DFTB is able to sufficiently improve the AMI 
potential energy surface. Whereas the frequency opti- 
mized variant of the SCC-DFTB method (SCC-DFTB- 
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MWM) improves the EIE in the HA, it does not retain 
the SCC-DFTB accuracy in the potential surface scans. 
Hence we decided to use the AMI and SCC-DFTB poten- 
tials for PIMD calculations. To conclude, AMI performs 
very well in HA, but it cannot properly describe poten- 
tial surfaces of the two rotational motions realized during 
simulations. On the other hand, SCC-DFTB gives worse 
results in HA, but it reproduces both barriers very well. 



FIG. 8: Relaxed potential energy scan of the vinyl group rota- 
tion in s-trans (3z7)-penta-l,3-diene (compound 1). For MP2, 
B98, DFTB, SCC-DFTB and SCC-DFTB-MWM only posi- 
tions and potential energies of minima and transition states 
are indicated. 



* Electronic address: jiri.vanicek@epfl.ch 

1 W. C. Alston, K. Haley, R. Kanski, C. J. Murray, and 
J. Pranata, J. Am. Chem. Soc. 118, 6562 (1996). 

2 A. Anbar, A. Jarzecki, and T. Spiro, Geochim. Cos- 
mochim. Acta 69, 825 (2005). 

3 E. Gawlita, V. E. Anderson, and P. Paneth, Eur. Biophys. 
J. 23, 353 (1994). 

4 P. S. Hill and E. A. Schauble, Geochim. Cosmochim. Acta 
72, 1939 (2008). 

5 D. A. Hrovat, J. H. Hammons, C. D. Stevenson, and W. T. 
Borden, J. Am. Chem. Soc. 119, 9523 (1997). 

6 K. E. Janak and G. Parkin, Organometallics 22, 4378 
(2003). 

7 K. Kolmodin, V. B. Luzhkov, and J. Aqvist, J. Am. Chem. 
Soc. 124, 10130 (2002). 

8 C. Munoz-Caro, A. Nino, J. Z. Davalos, E. Quintanilla, 
and J. L. Abboud, J. Phys. Chem. A 107, 6160 (2003). 

9 M. W. Ruszczycky and V. E. Anderson, J. Mol. Struct. 
THEOCHEM 895, 107 (2009). 

10 M. Saunders, M. Wolfsberg, F. A. L. Anet, and O. Kronja, 
J. Am. Chem. Soc. 129, 10276 (2007). 

11 L. M. Slaughter, P. T. Wolczanski, T. R. Klinckman, and 
T. R. Cundari, J. Am. Chem. Soc. 122, 7953 (2000). 

12 V. V. Smirnov, M. P. Lanci, and J. P. Roth, J. Phys. Chem. 
A 113, 1934 (2009). 

13 A. Zeller and T. Strassner, Organometallics 21, 4950 
(2002). 

14 E. L. Pollock and D. M. Ceperley, Phys. Rev. B 36, 8343 
(1987). 

15 M. Boninsegni and D. M. Ceperley, Phys. Rev. Lett. 74, 
2288 (1995). 

16 L. Torres, R. Gelabert, M. Moreno, and J. Lluch, J. Phys. 
Chem. A 104, 7898 (2000). 

17 T. Ishimoto, Y. Ishihara, H. Teramae, M. Baba, and U. Na- 
gashima, J. Chem. Phys. 128, 184309 (2008). 

18 R. D. Bardo and M. Wolfsberg, J. Phys. Chem. 80, 1068 
(1976). 



19 L. I. Kleinman and M. Wolfsberg, J. Chem. Phys. 59, 2043 
(1973). 

20 J. Bigeleisen, J. Am. Chem. Soc. 118, 3676 (1996). 

21 D. A. Knyazev, G. K. Semin, and A. V. Bochkarev, Poly- 
hedron 18, 2579 (1999). 

22 R. D. Cowan, The Theory of Atomic Structure and Spectra, 
Los Alamos Series in Basic and Applied Sciences (Univer- 
sity of California Press, 1981). 

23 W. V. Doering and X. Zhao, J. Am. Chem. Soc. 128, 9080 

(2006) . 

24 W. V. Doering and E. J. Keliher, J. Am. Chem. Soc. 129, 
2488 (2007). 

25 D. Chandler, Introduction to modern statistical mechanics 
(Oxford University Press, New York, 1987). 

26 R. Feynman and A. Hibbs, Quantum Mechanics and Path 
Integrals, International Series in Pure and Applied Physics 
(McGraw-Hill, 1965). 

27 C. Predescu, D. Sabo, J. D. Doll, and D. L. Freeman, J. 
Chem. Phys. 119, 10475 (2003). 

28 T. Yamamoto and W. H. Miller, J. Chem. Phys. 122, 
044106 (2005). 

29 J. Vamcek and W. H. Miller, J. Chem. Phys. 127, 114309 

(2007) . 

30 W. Wang and Y. Zhao, J. Chem. Phys. 130, 114708 (2009). 

31 J. Vamcek, W. H. Miller, J. F. Castillo, and F. J. Aoiz, J. 
Chem. Phys. 123, 054108 (2005). 

32 M. F. Herman, E. J. Bruskin, and B. J. Berne, J. Chem. 
Phys. 76, 5150 (1982). 

33 J. Vamcek and W. H. Miller, in Proceedings of the Eighth 
International Conference: Path Integrals from Quan- 
tum Information to Cosmology, edited by C. Burdik, 
O. Navratil, and S. Posta (JINR Dubna, 2005). 

34 C. Predescu, Phys. Rev. E 70, 066705 (2004). 

35 D. Case, T. Darden, I. T.E. Cheatham, C. Simmerling, 
R. D. J. Wang, R. Luo, K. Merz, D. Pearlman, M. Crowley, 
R. Walker, et al., Amber 9 (2006), University of California, 
San Francisco. 



14 



36 D. Case, T. Darden, I. T.E. Cheatham, C. Simmerling, 
J. Wang, R. Duke, R. Luo, M. Crowley, R.C.Walker, 
W. Zhang, et al., Amber 10 (2008), University of Cali- 
fornia, San Francisco. 

37 D. Chandler and P. G. Wolynes, J. Chem. Phys. 74, 4078 
(1981). 

38 W. R. Roth and J. Konig, J. Liebigs Ann. Chem. 699, 24 
(1966). 

39 A. D. Becke, J. Chem. Phys. 107, 8554 (1997). 

40 H. L. Schmider and A. D. Becke, J. Chem. Phys. 108, 9624 
(1998). 

41 M. J. S. Dewar, E. G. Zoebisch, E. F. Healy, and J. J. P. 
Stewart, J. Am. Chem. Soc. 107, 3902 (1985). 

42 M. Elstner, D. Porezag, G. Jungnickel, J. Eisner, 
M. Haugk, T. Frauenheim, S. Suhai, and G. Seifert, Phys. 
Rev. B 58, 7260 (1998). 

43 J. Merrick, D. Moran, and L. Radom, J. Phys. Chem. A 
111, 11683 (2007). 

44 H. A. Witek and K. Morokuma, J. Comput. Chem. 25, 
1858 (2004). 

45 E. Malolepsza, H. A. Witek, and K. Morokuma, Chem. 
Phys. Lett. 412, 237 (2005). 

46 T. Kriiger, M. Elstner, P. Schiffels, and T. Frauenheim, J. 
Chem. Phys. 122, 114110 (2005). 

47 H. Flyvbjerg and H. G. Petersen, J. Chem. Phys. 91, 461 
(1989). 

48 B. J. Berne and D. Thirumalai, Annu. Rev. Phys. Chem. 
37, 401 (1986). 

49 G. J. Martyna, M. L. Klein, and M. Tuckerman, J. Chem. 
Phys. 97, 2635 (1992). 

50 J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, and 
D. A. Case, J. Comput. Chem. 25, 1157 (2004). 

51 M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, 
M. A. Robb, J. R. Cheeseman, J. A. Montgomery, Jr., 
T. Vreven, K. N. Kudin, J. C. Burant, et al., Gaussian 03, 
Revision E.01, Gaussian, Inc., Wallingford, CT, 2004. 

52 B. Aradi, B. Hourahine, and T. Frauenheim, J. Phys. 
Chem. A 111, 5678 (2007). 

53 D. E. Sunko, K. Humski, R. Malojcic, and S. Borcic, J. 
Am. Chem. Soc. 92, 6534 (1970). 

54 J. C. Barborak, S. Chari, and P. v. R. Schleyer, J. Am. 
Chem. Soc. 93, 5275 (1971). 

55 J. J. Gajewski and N. D. Conrad, J. Am. Chem. Soc. 101, 
6693 (1979). 

56 T. Hascall, D. Rabinovich, V. J. Murphy, M. D. Beachy, 
R. A. Friesner, and G. Parkin, J. Am. Chem. Soc. 121, 
11402 (1999). 

57 B. R. Bender, J. Am. Chem. Soc. 117, 11239 (1995). 

58 F. Abu-Hasanayn, K. Krogh-Jespersen, and A. S. Gold- 
man, J. Am. Chem. Soc. 115, 8019 (1993). 

59 D. Rabinovich and G. Parkin, J. Am. Chem. Soc. 115, 353 
(1993). 

60 W. R. Angus, C. R. Bailey, J. B. Hale, C. K. Ingold, A. H. 
Leckie, C. G. Raisin, J. W. Thompson, and C. L. Wilson, 
J. Chem. Soc. 62, 971 (1935). 

61 O. Redlich, Z. Phys. Chem. Abt. B 28, 371 (1935). 

62 J. J. P. Stewart, J. Comput. Chem. 10, 209 (1989). 

63 G. B. Rocha, R. O. Freire, A. M. Simas, and J. J. P. Stew- 
art, J. Comput. Chem. 27, 1101 (2006). 

64 J. P. McNamara, A.-M. Muslim, H. Abdel-Aal, H. Wang, 
M. Mohr, I. H. Hillier, and R. A. Bryce, Chem. Phys. Lett. 
394, 429 (2004). 

65 M. P. Repasky, J. Chandrasekhar, and W. L. Jorgensen, 
J. Comput. Chem. 23, 1601 (2002). 



J. Stewart, J. Mol. Model. 13, 1173 (2007). 
J. J. P. Stewart, Mopac 2007 (2007), Stewart Computa- 
tional Chemistry, Colorado Springs, CO, USA. 
Y. P. Liu, G. C. Lynch, T. N. Truong, D. H. Lu, D. G. 
Truhlar, and B. C. Garret, J. Am. Chem. Soc. 115, 2408 
(1993). 

Namely, the barrier height of internal hindered rotation of 
the methyl group. 



